Iranian Journal of Numerical Analysis and Optimization 
Vol. 18, No. 1, 2023, pp 1-18 ao 


https: //doi.org/10.22067 /ijnao.2022.72900.1063 _ 

https://ijnao.um.ac.ir/ 

How to cite this article ro} 
BY NC Research Article 


Applying the meshless Fragile Points 
method to solve the two-dimensional 
linear Schrodinger equation on arbitrary 
domains 


D. Haghighi, S. Abbasbandy*and E. Shivanian 


Abstract 


The meshless Fragile Points method (FPM) is applied to find the numerical 
solutions of the Schrédinger equation on arbitrary domains. This method 
is based on Galerkin’s weak-form formulation, and the generalized finite 
difference method has been used to obtain the test and trial functions. For 
partitioning the problem domain into subdomains, Voronoi diagram has 
been applied. These functions are simple, local, and discontinuous poly- 
nomials. Because of the discontinuity of test and trial functions, FPM 
may be inconsistent. To deal with these inconsistencies, we use numerical 
flux corrections. Finally, numerical results are presented for some exam- 
ples of domains with different geometric shapes to demonstrate accuracy, 
reliability, and efficiency. 
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1 Introduction 


Numerical methods are mainly used to solve partial differential equations 
and have been studied, for example, the finite element method [1], finite vol- 
ume method [4], and boundary element method [13] to discretize the spatial 
dimension can be mentioned. In these methods, the accuracy of the method 
may be affected by deforming the elements or meshes. Therefore, meshless 
methods such as element free Galerkin [3] and meshless local Petrov-Galerkin 
[2] were considered. In these methods, the trial and test functions must be 
continuous, and usually, the trial functions in these methods are compli- 
cated. Dong et al. [6] introduced a new meshless method in which test and 
trial functions are considered as simple, local, and discontinuous polynomials. 
Very recently, this method has been used to solve the two-dimensional hy- 
perbolic telegraph equation [8]. This new method is called the Fragile points 
method (FPM), which we will study in this article. This method is also used 
for complex and irregular domains, which are discussed in this study. 


Solving the Schrédinger equation is very important in quantum dynamic 
calculations, and it has received a lot of attention as a model that describes 
several important chemical and physical phenomena [11]. This equation is 
derived from the vector wave equation for the electric field, which governs 
the propagation of electromagnetic waves in an inhomogeneous medium [10]. 
Schrédinger equations are also applicable in underwater acoustics, optics, 
and the design of optoelectronic devices [5]. 


We consider the two-dimensional time-dependent Schrédinger equation 
with the form 


Ou 5 
—t5 (6 2) = V~u(x, t) + w(x)u(x, t), x EQ, (1) 


with initial conditions 
u(x, 0) = g(x), (2) 
and the boundary conditions 
u(x,t) =hi(x,t), x€Tp, Vu.n(x, t) = ho(x,t), xETy. (3) 
In the above equation, w(x) is an arbitrary potential function. 


In the rest of this paper, in Section 2, the process of obtaining test and 
trial functions is described. In Section 3, the implementation of numerical 
flux corrections is given. Some numerical results and examples are provided 
in Section 4. Finally, the conclusions reached from using the FPM for the 
two-dimensional linear Schrédinger equation are expressed in Section 5. 
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2 Polynomial discontinuous trial and test functions 


In this section, we describe the process of obtaining local, simple, discontin- 
uous, polynomial, point-based trial and test functions. In order to divide the 
domain into subdomains, we consider several scattered points in the domain 
Q and its boundary 0. This subdivision should be such that each subdo- 
main contains only one point. Domain partitioning can be done in several 
ways, for example, the Voronoi diagram partition, quadrilateral and trian- 
gular partition (in 2D), tetrahedron and hexahedron partition (in 3D), and 
so on. In this study, the Voronoi diagram method has been selected. In the 
present FPM, only nonuniform or uniform points inside and on the domain 
boundary are applied, and it is a meshless method. 


The trial function up, in the subdomain Ep that includes the point Pp, 
can be written as 


un(x,t) = uo(x, t) + (k — Xo) Vu(x,t)|p,, x © Eo. (4) 


In the above equation, uo is the value of uz, at Po and xg denotes the coor- 
dinate of the point Pp. 


The gradient of Vu at Po is yet unknown. We employ the generalized 
finite difference method to calculate Vu at Po in terms of the values of up, 
at several neighboring points of Pp. We name these neighboring points as 
G1; 92;---,Gm- In the following, to calculate the amount of the gradient of 
Vu at Po, we minimize a weighted discrete L? norm J so that 


J= s ( Vulpy- (xi — Xo)” — (us —)) ws, (5) 
1=0 


where w; denotes the value of weight function at q;, x; is the coordinate 
vector of qg;, and u; is the value of up, at gq; (¢ = 1,2,...,m). For convenience, 
we assume that w is constant. Due to the stationarity of J, we have 


= 
Vu= (ATA) AY (Um — Uol mn) , (6) 
where 
L1—Xo Yi — Yo Uy 1 
T2—2X9Q Y2— Yo U2 1 
= ’ Um = , Ln = 
Tm — LQ Ym — Yo Um 1 mx1 


Also equation (6) can be expressed at point Pp as follows: 


Vu = Bug, (7) 
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where 
—-110...0 uo 
= — , U1 
B = (A7A) ae : up = 
—-10--- 01 Um 


mx(m+1) 


Also by substituting (7) into (4), the relation between up, and ug will be 
obtained as 


up =Nug, forallxe€ Ey, N= [x—xo]B+[1,0,...,0]1x(m4i)- (8) 


3 Numerical flux corrections 


We can rewrite Schrédinger equation (1)—(3) using mixed form as follows: 


o(x,t) = Vu(x, t), in Q, 
Ou 
—V.o(x,t) = ia, 2) + w(x)u(x, t), in Q, (9) 
u(x,t) = hi(x,t), in Tp, 
o.n(x,t) = ho(x,t), in Tn. 


By multiplying the first and second equations in (9) by the test functions + 
and v, respectively, and integrating it on the subdomain E, 


| o7;,°TAQ = ‘ Vun(x, t) -7dQ, (10) 

E E 

| -V -ovdQ = if Fel, Qua + | w(x)u(x, t)vdQ, (11) 
E gE Ot E 


using the Green formula and by summing these equations over all subdo- 
mains, we have 


J enrao = -{ unV.TdQ + ys 
Q Q 


| ay,n.rdD, (12) 
Beg” IE 


| o7,.VvdQ. = S- 
Q 


| pnw +i Sp sted+ f w(xulx,t)vad. 
fen OE Q Ot Q 


(13) 
In the above equations, values 6, and ti, represent approximations op, and 


up, on OF. These values are named numerical fluxes. To simplify (12) and 
(13), we define the operators average and jump, where by these operators, 
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we can manage the numerical fluxes. As regards, T =T,+Tpn+Tw (Lp is 


the set of all internal boundaries of subdomains), using [6, Table 3.1 ], and 
substituting the interior penalty numerical fluxes, we have 


ee S- [van [v] + {Vv} [un]) a 


EEQ e€l,Ulp ~ © 
+ b> 7 fw [un] dl = S- [ (e- ven) hi (x, t)aT 
e€l,Urp e€l'p 


+ i: w(x)u(x,t)vdX+ S° | vho(x,t)dP 


Q ecIn ° © 
+i f Sextuan. 


The above equation is the formula of FPM, which is called FPM-primal 
method. 


The method (matrix form) can be expressed as follows: 
(K —W)u-iCu=F. (14) 
Using 0-weighted scheme [12], the above equation can be written as follows: 


yett _ u” 
(K — W)(6u"*? + (1 — @)u”) — iC( ie =r", (15) 


In (15), u*(x) = u(x, kAt), where At is the time step and 0 < 6 < 1. By 
substituting values B instead of Vv and Vu, N instead of uy; and vy in the 
formula of FPM-primal, the point stiffness matrices K, C, W and also the 
right vector F will be achieved as follows: 


w= / N’Nw(x)dQ, c= | N7NaQ, Kp= | B’BdQ, E£€Q, 
E E E 


(16) 
=f 
K,= > [otntn, +N? n,B,)dP + 7. [ NtNiar 

=i 
+ | (BInNo +NjnsBy)dP + 7 ; N2NodP 

a4 
+ | (BEnTN, +NjnsB,)d0 + - / N? Nod? 

= 
+— | (Bin}N2 + N/m Bo) + 2 [ NtNuar, e € OE, NOE2, 
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Kp = - [e'ntn +N’nB)dr + . [NtNar, e€Tp, (17) 


and it can also be written 


Fy = [NT habe oar, e€Ty, 
Fp= [gay ~B'n)hi(x, t)dD, e€Tp. 


4 Numerical results 


In this section, we will study some numerical examples, and using the results 
obtained from the application of FPM on these examples, the accuracy and 
stability of the method are investigated. All examples were done in MATLAB 
software on a Core i5, 2.67 GHz CPU machine with 4 GB of memory. The 
relative errors used in this section are defined as follows: 

7 Jun — ull z2 _ |Vupn — Vull 2 


To = 1: 
ullc2 |Vullz2 


Also we calculate the convergence orders in space and time via 
log 10 (2) log 10 (=) 
7 , C—order(time) = —— 
log 10 (=) log 10 (=) 
12 


such that h, and At; correspond to e; and also hg and At, correspond to 
error €2. In numerical examples, we consider errors e; and € as relative error 
TO. 


C — order(space) = 


Example 1. (a) We first consider (1) with potential function w(x, y) = 1— 


2p and exact solution u(a, y,t) = e’'x?y? in the region Q = (0, 1] x [0, 1]. 
Boundary and initial conditions are obtained using the exact solution [7, 11]. 
In Table 1 relative errors are shown for the number of different points of 
the domain that are uniformly distributed and all boundary conditions are 
considered Dirichlet. This table shows the good accuracy and stability of the 
method, and as the number of points increases, the accuracy of the method 
improves. In addition, as shown in Figure 1, the relative errors decrease as 
the number of points increases. 
(b) Next, we consider the boundary conditions as follows: 


u(0,y,t) = u(#,0,t) =0, Vu.n(1,y,t) = 2ey?, u(x,1,t) = e*x?. 
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Table 1: Relative errors of the method for example l(a) at T = 1 and At = 0.05 with 
6 = 0.6 and uniform points. 


h Number of parameters Errors CPU time (s) C-order 
point 

0.25 N =25 he = 0.1 3.942832 x 10~? 0.31 - 
n= 2.5 2.523913 x 107+ 

0.1 N=121 he = 0.1 5.73321 x 10-3 0.63 2.1044 
n=4 3.52258 x 107? 

0.04 N=676 he = 0.1 9.91802 x 10-4 9.34 1.9148 
n=9 9.26574 x 10-3 

0.02 N = 2601 he = 0.1 4.30894 x 10-4 163s 1.2027 
n=9 4.27138 x 1078 


In Table 2, relative errors have been reported for the number of different 
points that are uniformly and nonuniformly distributed over the domain. 
Comparing the results, we find that in this example, how the points are dis- 
tributed does not have much effect on the accuracy of the method. Also 
Figure 2 shows error plots for N = 676 uniform and nonuniform points. Fig- 
ures 3 and 4 demonstrate the appropriate accuracy of the FPM for different 
times for x = 0.6, 0 = 0.51, At = 0.01, N = 121,h, = 0.1, and n = 4. 
Therefore, according to these results, the method has appropriate accuracy 
for different boundary conditions and is also convergent. In addition, as the 
final time increases, accuracy is maintained and FPM is stable. 

Compared to [15] and [16], the proposed method achieves almost the same 
accuracy in much less time. 


log 0(r) 


-3.5 


log1 O(N) 


Figure 1: Relative errors for Example 1 at T = 1 and At = 0.05. 
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Table 2: Relative errors of the method for Example 1(b) at T = 1, 6 = 0.51 and At = 0.01 
for points with uniform and nonuniform distribution. 


uniform points 


h Number of points Errors CPU time (s) C-order 

0.25 N =25 6.476626 x 10-? 0.40 : 
2.636926 x 107! 

0.1 N=121 5.763461 x 1073 1.08 2.1864 
3.419285 x 1072 

0.04 N = 676 2.552031 x 1073 27.9 0.8891 


1.216389 x 10~? 


nonuniform points 


Number of points Errors CPU time (s) 
N = 25 4.273177 x 107-2 0.33 
1.484748 x 107! 
N=121 9.238729 x 10-8 0.94 
9.161839 x 1072 
N = 676 2.423225 x 1078 o7.7 


1.216389 x 10~? 
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Absolute error Absolute error 


Figure 2: Comparison the absolute errors related to Example 1(b) for T = 1, At = 0.01, 
6 = 0.51, and N = 676 for uniform and nonuniform points. 


Imaginary part 


Exact solution 
0.3 ° FPM solution 


Figure 3: Comparison the imaginary parts of exact and numerical solutions related to 
Example 1 for z = 0.6, At = 0.01, and N = 121. 


Example 2. In this example, we consider (1) with N = 676 uniform points 
in the domain 2 = [0,1] x [0,1] such that 


Ag? + dy? — 42 — 4y + 6? — 4842 

ee 

(x—0.5)? (y—0.5)? 
B B 


w(z,y) = 


u(a,y,t) = exp( it). 
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Real part 


Exact solution 


o.2b fo) FPM solution 


L L L L L L L L 
oO 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 


Figure 4: Comparison the real parts of exact and numerical solutions related to Example 
1 for « = 0.6, At = 0.01 and N = 121. 


Table 3: Relative errors of the method for Example 2 with he = 0.1, 0 = 0.52, 7 = 11, 
and At = 0.01. 


Final time ro ry CPU time (s) 
S11 1.85272 x 10-+ 1.18535 x 107? 28.3 
T=5 3.84086 x 10-4 2.00094 x 1072 99.5 
T= 10 5.48035 x 10-4 2.10124 x 107? 201.3 
T=15 5.00376 x 10-4 2.07114 x 107? 281.5 


With Dirichlet boundary conditions for At = 0.01, 0 = 0.52, and 6 = 2, 
relative errors related to different final times are shown in Table 3. This table 
shows the stability of the method over time. In the following, Figures 5 and 
6 show the plots of imaginary and real parts of numerical and exact solutions 
for N = 2601 uniform points with h, = 0.1 and 7 = 11. These figures indicate 
that the method is also accurate for a large number of points. Also, the plot 
of errors for N = 676 uniform and nonuniform points is provided in Figure 
7. As this figure shows, under similar conditions, the error of the proposed 
method is less for points with a uniform distribution. 


Example 3. a) Now we solve the previous example in an L-shaped domain 
that has the following boundaries: 


Q, = {0} x [0,1], Q» = {0.36} x (0.52, 1], 
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Imag(Numerical) Imag(Exact) 


Figure 5: Comparison of imaginary parts of the numerical and exact solutions for Ex- 
ample 2 with 7 = 11, he = 0.1, = 0.52, At = 0.01, N = 2601, and T = 1. 


Real(Numerical) Real(Exact) 


Figure 6: Comparison of real parts of the numerical and exact solutions for Example 2 
with 7 = 11, he = 0.1, 9 = 0.52, At = 0.01, N = 2601, and T = 1. 


Q3 = [0,1] x {0}, 4 = [0.36, 1] x {0.52}. 


Figure 8 shows the uniform distribution of N = 484 points in this domain. If 
we consider the boundary conditions completely Dirichlet, then we have Table 
4 for the number of different points of the domain that are uniformly selected. 
As this table shows, according to the number of points, the numerical results 
have good accuracy that is obtained in a short time. Figure 9 shows the 
relation between the distance between points with uniform distributions with 
relative errors. 


b) Next, we consider a circular domain as follows: 


Q= {(e,y) ER: Ver +y <1}, 
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Absolute error Absolute error 


x10 #408 


Figure 7: Plot of error for Example 2 based on uniform points for 0 = 0.52 N = 676, 
T =1, and At = 0.01. 


and the equation is solved by FPM. Table 5 shows the relative error values 
and convergence orders over time. According to this table, there is no need to 
reduce the time step too much, because making it smaller does not have much 
effect on accuracy, and we should improve the accuracy by changing other 
parameters or the number of points. Also, due to the circular amplitude and 
the number of points used, the relative errors obtained are acceptable and 
are obtained in a short time. 


. 
' 
oe 
' 
‘ 


tot 
eee eee 
eee ee nenee 


Pe eR 
1 
’ 


+ }+++4+1++1 


Figure 8: L-shaped domain related to Example 3(a) with N = 484 uniform point: 


e 


Example 4. In this example, we consider the time-dependent Schrédinger 
(1)-(3) in (a,y) € Q = [0,z] x [0,7] with initial condition u(z,y,0) = 
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Table 4: Relative errors of the method for Example 3(a) at T = 1 and At = 0.01, 
0 = 0.65. 


Number of points parameters Errors CPU time(s) 
N =49 he=1 ro = 4.796603 x 10-3 0.69 
n= 75 r, = 1.608921 x 107! 
N=121 he =0.001 rp = 3.516616 x 10-3 1.20 
7 = 190 r, = 1.044712 x 1071 
N = 484 he =0.1 ro = 3.649402 x 10-4 12.60 
n= 75 ry = 2.631145 x 10-? 
N = 1849 he =0.01 rp = 3.196693 x 10-4 280.72 
n = 27 r= alias «10-7 
-0.5 
-—1PF Sa 
ash 
-ash 
=3-+ 
S87 16° 15 —1.4 4 1 -0.9 -0.8 


1.3 -1.2 
log10(h) 


Figure 9: Error curves with respect to the distance of selected points from the domain 
to each other for example 3(a) with At = 0.01 and T = 1s. 


Table 5: Relative error values and convergence orders over time for example 3(b) with 
N = 529 uniform points, T = 1, 6 = 0.65, he = 0.01, and 7 = 16. 


Time step To ry CPU time(s) C-order 
At=0.08 1.261254 10-2 6.019296 x 10-2 6.88 - 

At =0.04 7.115889 x 10-3 4.330438 x 107? 9.04 0.8257 
At =0.02 3.882812 10-3 3.440278 x 10-2 13.67 0.8739 
At=0.01 3.042351 x 107% 3.105425 x 107? 23.26 0.3519 


sin(x) sin(y) and Dirichlet boundary conditions that are zero on all sides, 
with the given potential function as w(x, y) = 3,(a,y) € Q. The analytical 
solution is as u(x, y,t) = e” sin(x) sin(y). 
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As you can see in Table 6, for a number of different points, the method 
has the appropriate accuracy to solve this example. Also, Figures 10 and 11 
demonstrate the plots for he = 0,1, 7 = 3.8, and N = 676 nonuniform points 
for T = 1s and T = 3s, respectively. These figures show the accuracy of 
FPM for the case where the points are considered nonuniform. Compared to 
[9], the proposed method reports better computational times and accuracy. 


Table 6: Relative errors of the method for Example 4 at T = 1, 6 = 0.52, and At = 0.01. 


h Number of parameters Errors CPU time (s) C-order 
points 

0.25 N = 25 he =0.1 1.797325 x 10-? 0.41 - 
n=1 1.185351 x 107! 

0.1 N = 121 he =0.1 1.89251 x 10-8 0.97 2.4566 
n= 2.7 1.92942 x 10~? 

0.04 N =676 fe=0.1 6.002710 10-7 9.47 1.2532 
n=3.8 4.150509 x 1073 

0.02 N= 2601 fe=O0.1 4.561287 x 10-7 107 0.3962 


n=2 6.729947 x 1073 


Real(Numerical) Real(Exact) 


Figure 10: Plots for real parts of the numerical and exact solutions related to Example 
4 for dt = 0.01, N = 676, T = 1, 6 = 0.54, 7 = 38, and he = 1 


Example 5. For the last example, we consider Schrédinger equation with 
the following exact solution and initial conditions: 


u(x, y,t) = e—) (sin(x) + cos(y)), u(x, y,0) = (sin(x) + cos(y)) . 


This equation is solved using the proposed method on the connected amoeba- 
like domain according to Figure 12 with N = 100 points that are nonuni- 
formly distributed in the domain. Table 7 shows relative errors and CPU 
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Imag(Numerical) 


Imag(Exact) 


Figure 11: Plots for imaginary parts of the numerical and exact solutions related to 
Example 4 for At = 0.01, N = 676, T = 5s, 0 = 0.54, 7 = 38, and he = 1. 


time related to different final times. Due to the nonuniform points and the 
domain of the problem in this example, the accuracy of the method is ac- 
ceptable. Also, plots related to numerical and exact solutions for T = 2s 
are presented in Figure 13, which confirms the suitability of the method for 
irregular domains. 


oa 


ee ae . 
A EN eciahe ae 
oO L L L ss L 1 
oO o.1 0.2 0.3 O.4 0.5 0.6 o.7 0.8 0.9 1 
x 


Figure 12: Domain of the problem in Example 5 with N = 100 selected points that are 
nonuniformly distributed. 


Table 7: Relative errors of the method for Example 5 for T = 1,5, 10,15, 20, At = 0.0097, 
and 6=0.5. 


T Parameters ro ry CPU time (s) 
1 he =0.001,7=215 5.697623 x 10-% 9.610123 x 107? 1.06 
5 he =0.001,7=600 6.276758 x 107° 9.613064 x 107? 3.22 
10 he =0.001 ,7=500 6.286391 x 10- 9.620433 x 10-2 6.05 
15 he =0.001,7=550 6.282664 x 10- 9.618269 x 10-2 8.69 
20 he =0.001 , 7 =600 6.279557 x 107° 9.616465 x 107? 11.33 
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Real(Numerical) Real(Exact) 


Imag(Numerical) Imag(Exact) 


Figure 13: Plots related to numerical and exact solutions for Example 5 with 6 = 0.5, 
he = 0.001, 7 = 50,T = 2s and At = 0.0097. 


5 Conclusion 


In this paper, the meshless Fragile Points Method (FPM) is used to obtain 
numerical solutions to the two-dimensional linear Schrédinger equation. This 
method is based on Galerkin’s weak form, and the test and trial functions 
are considered simple, local, and discontinuous polynomials. Numerical flux 
corrections have been used to deal with inconsistencies due to the discon- 
tinuity of trial functions. Finally, the efficiency, stability, and accuracy of 
the method were evaluated with several numerical examples. In these nu- 
merical examples, the accuracy and stability of the method were evaluated 
both for a large number of points and for the case where the points were 
selected nonuniformly. We also got good solutions for larger final times and 
the problems with irregular domains. 


According to the results of the tables and comparison of the curves ob- 
tained by FPM with the exact curves, it can be seen that the method is stable 
and has good accuracy. Also, the method does not have much computational 
cost and depending on the number of points used, it will achieve numerical 
solutions with good accuracy in a short time, which is an advantage over 
the finite element. Other advantages of this method over other numerical 
methods are described in detail in [14, Table 1 ]. 
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